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ABSTRACT 

The evolution of globular clusters due to 2-body relaxation results in an outward flow of en- 
ergy and at some stage all clusters need a central energy source to sustain their evolution. 
Henon provided the insight that we do not need to know the details of the energy production 
in order to understand the relaxation-driven evolution of the cluster, at least outside the core. 
He provided two self-similar solutions for the evolution of clusters based on the view that 
the cluster as a whole determines the amount of energy that is produced in the core: steady 
expansion for isolated clusters, and homologous contraction for clusters evaporating in a tidal 
field. The amount of expansion or evaporation per relaxation time-scale is set by the instan- 
taneous radius or number of stars, respectively. We combine these two approximate models 
and propose a pair of Unified Equations of Evolution (UEE) for the life cycle of initially 
compact clusters in a tidal field. The half-mass radius increases during the first part (roughly 
half) of the evolution, and decreases in the second half; while the escape rate approaches a 
constant value set by the tidal field. We refer to these phases as 'expansion dominated' and 
'evaporation dominated'. These simple analytical solutions of the UEE immediately allow us 
to construct evolutionary tracks and isochrones in terms of cluster half-mass density, cluster 
mass and galacto-centric radius. From a comparison to the Milky Way globular clusters we 
find that roughly one-third of them are in the second, evaporation-dominated phase and for 
these clusters the density inside the half-mass radius varies with the galactocentric distance 
i?G as ph C)c Rq'- The remaining two-thirds are still in the first, expansion-dominated phase 
and their isochrones follow the environment-independent scaling ph c< A'P, where M is the 
cluster mass; that is, a constant relaxation time-scale. We find substantial agreement between 
Milky Way globular cluster parameters and the isochrones, which suggests that there is, as 
Henon suggested, a balance between the flow of energy and the central energy production for 
almost all globular clusters. 
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1 INTRODUCTION 

Globular clusters strive their entire life for thermal equilibrium, but 
their negative heat capacity prevents them from ever reaching this. 
As a result there is a continuous flow of energy that is conducted 
outwards by relaxation. Any equations attempting to capture clus- 
ter evolution could be very complex and non-linear with many vari- 
ables to keep track of since the dynamical evolution of star clusters 
is the result of several processes, including two-body relaxation, in- 
teractions with binary stars, escape across the tidal boundary, and 
the internal evolution and mass-loss of single and binary stars. Fig[T] 
illustrates typical results, in terms of the evolution of the core ra- 
dius, the half-mass radius and the tidal radius. Within the first few 
times 10 Myr there is a rapid phase of mass segregation, leading 
to contraction of the core. (Note the logarithmic time scale in the 
figure.) At the end of this phase the half-mass radius starts to in- 
crease, and the expansion continues, but decelerates, over the next 



8 Gyr approximately. At the same time the cluster is losing mass, 
both by the escape of stars and by the effects of stellar evolution, 
and this causes a decrease in the tidal radius. For the remaining evo- 
lution (up to 12 Gyr) the half-mass radius also contracts. Except for 
the initial phase of mass segregation, the core radius remains small 
compared with the half-mass radius; it also begins its contraction at 
about 6 Gyr, i.e. earlier than the half-mass radius, and subsequently 
drops to and fluctuates around a lower value. 

Our theoretical aim in this paper is to provide a physically mo- 
tivated and simple prescription for the behaviour of the half-mass 
and tidal radii in the period following the end of the phase of mass 
segregation, i.e. all except about the first 1% of the evolution of the 
cluster. With this we are able to construct evolutionary tracks and 
isochrones for clusters evolving in a tidal field and provide a theo- 
retical framework for empirically established correlations between 
structural parameters and t heir environment as found for Milky 
Way globular clusters (e.g. iDiorgovskil 1 19951 ; iMcLaughliiJ |200(T) 
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Figure 1. The evolution of the tidal, half-mass and core-radii (from top to 
bottom) in the evolution of a Monte Carlo model of the globular cluster M4 
(from .Heggie & Giersz„2008i) . 



and f o r extra-galactic globular cluster systems (e.g . Ijordan et"al 
20051: iBarmbv et ai] |2007| ; iMcLaughlin et al.l |2008| ; iHarris et al 
201(1) . We do not aim to explain the shape and dependence on envi- 



ronment of the globular cluster mass function at this stage, however. 
The other principal exclusion is that we do not attempt to model the 
evolution of the core. The fact that this programme will succeed is 
due to a number of important discove ries by Henon . 

The first of these discoveries jHenonl Il96 ill was a simple 
model of an idealised star cluster evolving in a tidal field. The 
main idealisations were that all stars have the same mass and do 
not evolve internally, the tidal field is modelled as a cutoff at the 
tidal radius, the system is non-rotating and spherically symmetric, 
there is no significant binary population, the velocity distribution 
of the stars at a point is isotropic, and gravitational encounters are 
modelled as a Fokker-Planck equation. Then Henon showed that 
the system could evolve at constant mean density, with stars escap- 
ing across the tidal boundary, and the structure of the cluster at any 
time is simply a scaled version of the structure at any other time, 
i.e. the evolution is 'self-similar'. In particular, this means that the 
half-mass radius is a constant fraction of the tidal radius. We shall 
adapt this model to describe the evolution at late times in models 
such as that shown in Fig[T] 

Henon's next advance in this field ( Henon|[l965il was the dis- 
covery of a similar model for an isolated system, i.e. one with no 
tide. It is again self-similar, but the total mass remains constant, and 
the half-mass radius increases. This we shall adapt to describe the 
evolution of models such as that shown in Fig[T] from the end of 
the phase of mass segregation until the point where the tidal radius 
begins to have a significant effect on the expansion of the half-mass 
radius. A further part of our adaptation will be to construct a way 
of bridging the transition between the two models, and to do so in 
a physically well motivated way. 

The third discovery of Henon explains why this programme 
can work even without detailed modelling of the core. In both of 
Henon's models the energy of the system increases, and Henon 
himself considered that the formation and hardening of binaries in 



the core was a plausible source of this energjQ. It might have been 
thought, therefore, that the evolution of the system would depend 
on details of this source of energy and the radius and density of the 
core. But a correct understanding of the connection between the 
core an d the evolution of the system as a whole came from the in- 
sight o f lHenonl ( fT97l . He realised that the rate of flow of energy is 
controlled by the system as a whole, and not by the core. That is, it 
is not the cluster that responds to whatever happens in the core, but 
the other way around. 

In Henon's picture the mechanism of energy generation is self- 
regulatory, as is the case in stellar interiors, where the luminosity 
of the star is determined by how much radiation can be transported 
through the envelope. The nuclear fusion rate in the core of the 
star adjusts to it. The application of this idea to stellar dynam- 
ics was a breakthrough allowing modellers to overcome the core 
collapse phase. In the Monte Carlo models of i HenoiJ ( 1 19751) the 
energy source is entirely artificial and unspecified. Later, it was 
found from A'^-body simulations of the long term (post-collapse) 
evolution of equal mass clusters that hard binaries act as the en- 
ergy source duiersz & Heggid 1 1994 iBaumgardt et"^ l2002h . In 
more realistic A'^-body simulations of star clusters these binaries 
are usually considered to be primordial, but other mechanisms of 
energy generation have been modell ed, including the action of a 
central intermediate-mass black hole dBavungardt et al.|[2004h and, 
less co ntroversially, mass-loss from stellar evolution. Gieles et^ 
( l2010h showed that mass loss as a consequence of stellar evolu- 
tion also leads to a self-regulatory energy production, which works 
together with hard binaries in driving the long term evolution. In- 
deed, it may even domi nate, as in the specific example of 47 Tuc 
dGiersz & He^l201lh . a high-concentration cluster in which the 
evolution of the core- and half-mass radii appears to be little af- 
fected by the primordial binary population. 

These self-regulatory mechanisms of energy generation take 
time to establish the balance between the energy generated in the 
core and the energy requirements of the overall evolution of the 
cluster. In Fig[T] this balance is reached at about the end of the 
phase of mass segregation. Mass segregation ceases close to the 
time when the core radius reaches a minimum for the first time. 
This was found from numerical simulations of clusters with a mod- 
erately wide mass spectrum (Giersz & Heggie 1996) a nd for clus- 
ters with a full mass spectrum (iPortegies Zwart et alj|2007h . The 
reason for this is n ot that equipartition has been reached, due to 
Spitzer's instability l lSpitze r 1969). The origin of this quasi steady- 
state distribution of stars of different masses is not well understood 
theoretically, but it is taken as a given here. We refer to the subse- 
quent evolution as 'balanced'. In much research on cluster dynam- 
ics this kind of evolution is usually associated with 'post-collapse' 
evolution, but this term is ambiguous; in Fig[T] for instance, the 
phrase 'post-collapse' might be used for the entire evolution after 
the end of mass segregation, but others might apply it to the evo- 
lution following the late decrease in the core radius in the last few 
Gyr. For this reason we prefer the term 'balanced evolution'. 

iGieles etal] ( |20I(]|) showed that the evolution of most glob- 
ular clusters is balanced. Their result was based, however, on a 
comparison of the behaviour of clusters expanding in isolation 
to the structural parameters of relatively massive globular clus- 
ters (>fewx 10'' Mq) and ultra-compact dwarf galaxies. It was as- 
sumed implicitly that these are little affected by the tide. In the 



' In fact lHenor] )l96ll) introduces the concept of energy generation by the 
core in the other way around: the core absorbs negative energy. 
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present study we add the effect of a tidal field, which slows down 
the expansion of all clusters at some point in their evolution. 

In Section[2]we discuss how the behaviour of the Henon mod- 
els in isolation and in tidal fields can be related to the flow of en- 
ergy, and we show how the two extremes, expansion without mass 
loss and homologous contraction in a tidal field, may be unified. 
Isochrones in physical units at an age of a Hubble time are extracted 
and are compared to the data of Milky Way globular clusters in Sec- 
tion[3] A summary and discussion are presented in Section|4] The 
specific symbols used in this study are summarised in Table [T] and 
details of the derivations are presented in Appendices. 



2 A UNIFIED MODEL OF EXPANSION AND 
EVAPORATION 

2.1 The flux of energy 

Our attempt to unify the evolution of mass and radius of the two 
models of Henon begins with one of the physical properties which 
the two models have in common, i.e. a flux of energy at the half- 
mass radius which is fed by an energy source in the core. We con- 
strain ourselves to the energy flow at this radius, because we can 
then construct a relatively simple set of relations for the behaviour 
of the bulk properties of the cluster. A full unification of the two 
models of Hen on req uires numerically solving the Fokker-Planck 
equations (e.g. iCohnlll979 . and follow-up studies). This would be 
needed to get the distribution function of energies and the density 
profile as a function of time, which is not the aim of this paper. For 
the present we adopt Henon's idealisation of systems in which all 
stars have the same mass m. Then we estimate the energy as usual 
by 



E = ~a- 



(1) 



where A'^ is the number of stars, rh is the half-mass r adius, and a is 
a 'form factor' for which the value 0.2 is often taken i lSpitzeJl987l. 
p. 12). 

In Henon's isolated model l lHenoij|l965h . all quantities on the 
right-hand side of equation ^ are constant except rh, and it follows 
from Henon's equations (3), (6) and (7) that 



E 

W\ 



rh 



0.0926 



(2) 



except that one also requires the value for the dimensionless half- 
mass radius R, which is given in Henon's paper on p. 64. Also, we 
have expressed the result in terms of the conventional half-mass 
relaxation time jSpitzeJ 19871) . i.e. 



Trh 



0.138- 



JVl/2 



.3/2 



(3) 



where rfi is the mean stellar mass (i.e. m in this case) and In A is the 
Coulomb logarithm, which we have equated with the factor Log n 

in Henon's equation (7). 

For Henon's tidally limited model ( lHenon|[l96 ll) . both A'^ and 
Th vary because the cluster loses stars at a constant density, and a 
similar calculation leads to the result that 



E 



0.0743 

Trh 



(4) 



What is noticeable about equations ^2} and © is how similar 
the final numerical coefficients are. Indeed, in constructing a uni- 
fied approximate model which includes the transition from nearly 



isolated evolution to tidally limited evolution, we shall make the 
assumption that the numerical coefficient is constant, i.e. that 



E 

W\ 



Trh ' 



(5) 



where ( ~ 0.08 in the case of equal masses. The accuracy of this 
approximation is comparable with the common approximation of 
treating a in equation l[T) as constant. 

Before proceeding further, we shall change one of the vari- 
ables in which the total energy E is expressed, because this will fa- 
cilitate one particular step in the further development of our model. 
Instead of using the half-mass radius, rh, we shall use the (half- 
mass) crossing time of the system, defined here as 

Tcr = (Gph)-'/' , (6) 

with ph = 3i\f/(87rrh), the cluster density within the half-mass 
radius, and AI = mN the total cluster mass. (Note that this differs 
from the conventional crossing time by a numerical factor; see also 
Table[T]) Putting this together with our assumption about the energy 
flux, we replace equations ^ and Q by0 



E 

Wl 



A. 

Trh 



5 iV 2 fcr 
"3iV 3~r 



(7) 



The constants become obvious when we look at the expression for 
the total energy (equation [TJ in terms of A'^ and Tcr (equation |6](: 
\E\ oc {mN)^^^ /tc/^ ■ In the next subsection the relative contribu- 
tions of expansion and evaporation to the energy flow are specified. 

2.2 The relative importance of expansion and evaporation 

Now we have a functional form that relates evaporation and expan- 
sion to the flow of energy (equation [7} we need to specify what 
the relative importance is of these two processes. The equations 
that follow (equaions [8[ll2t form the base for our analysis. We 
introduce two new coefficients, ^ and Xi to quantify the (dimen- 
sionless) evaporation rate and the (dimensionless) expansion rate, 
respectively. 



and 



N ' 



(8) 



(9) 



(10) 



so that (see equation|7} 

. 2 5, 

For Henon's isolated model, ^ = 0, while x — for the tidally 
bound model. 

To get a full description of the behaviour of the mass and 
half-mass radius of the cluster we need to relate one of the co- 
efficients ^ or X to Tcr and A'^. Here we focus on the dimension- 
less escape-rate 5. By estimating the fraction of stars above the es- 
cape velocity for tidally limited clusters with different ratios rh/rj, 
wher e r j is the Jacob i (tidal) radius of the cluster, iLeei iim^) and 
iGieles & Baumgardj (12008) have shown that In^ oc r^/r], ap- 
proximately. In principle we could proceed with this exponential 
function for ^, but a more convenient approximation will make 

^ For clusters with a stellar mass spectrum there will be an additional term 
— {5/3)(7Ti/?fi) on the right-hand side. 
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Table 1. Overview of the specific symbols used in this study. The symbols with the subscript 1 are used to refer to the maximum value the variable can have. 



Symbol Definition Description 



Ph 


3M/{8TTrl) 


Mean density within the half-mass radius 


P.I 


3M/(47rr3) 


Mean density within the Jacobi radius 






Maximum ratio of the half-mass radius over the Jacobi radius 






A measure of the crossing time at the half-mass radius 


r-I 

cr 


(Gpj)-i/2 


Crossing time at the Jacobi radius 


Tcrl 


^/2K/rJ)?/VJ 


Crossing time at the half-mass radius for a homologous cluster, equation jl6) 


Tih 


CNTcr 


Approximate half-mass relaxation time, with C ~ 0.006 (equations 1 1 3 l&[T4l 


TevO 


see equation i23\ 


Total life-time of a cluster of A^o stars 




see equation (23) 


Remaining Hfe-time, or 'life-expectancy', of a clusters of stars 


c 




Fraction of the total energy conducted per r^ij 


X 




Dimensionless expansion rate 




-(7V/Ar)r,h 


Dimensionless escape rate 



what is to follow much easier. In fact, iGieles & Baumgardtl l20o3) 
show that ^ oc (rh/rj)^/^ provides a satisfactory approximation to 
the exponential function (in the relevant range of rh/rj). This rela- 
tion can be expressed equivalently as f oc Ta. As a cluster expands, 
Tcr increases, until the point where the cluster evolution becomes 
tidally limited, i.e. t he point where its evolution begins to be analo- 
gous to the model o f lHenonl(ll96ll) . Let Tcrl denote the value of Tcr 
during this part of the evolution. Since Tcv is now constant, x ~ ^ 
(by equation|9]l, and so equation l IlQt implies that ^ = {3/5)( when 
Tcr = Tcrl . Putting these results together with equation JlOb . we ar- 
rive at a pair of Unified Equations of Evolution (UEE) 



X 



Nt, 



rh 



N 

'7'cr "^rh 3 



3 ^ Tcr 
5 Tcrl 



Tcr 
Tcrl 



(11) 

(12) 



which we adopt henceforth. These two expressions for ^ and x 
capture the behaviour of clusters for all values ^ Tcr ^ Tcri. 
For Tcr ^ Tcrl they give the desired self-similar expansion (x — 
(3/2)(^ = constant and ^ ~ 0) and the growth of Tcr stops when 
Tcr = Tcrl (x " and 5 = (3/5)0. Iri this regime we assume that 
the radial scale of the cluster is set by the tidal radius, and so the 
ratio of the half-mass to tidal radii takes a constant value, w hich we 
denote by (rh/rj)i. In his equal-mass homological model IllenonI 
( Il96ll) found the value (rh/rj)i ~ 0.145. 

In Fig. [2] we show the behaviour of ^ and x in units of C, as 
a function of Tct/tcti- The top axis shows how rh/rj relates to 
Tcr/Tcri, wMch wlU be explained in more detail in the next subsec- 
tion (eguationllSt. 



2.3 Time scales 

Before proceeding to solve the equations of our model for clus- 
ter evolution, it is convenient to collect one or two straightforward 
formulae for future reference. First, we often approximate the half- 
mass relaxation time (equation[3]( by 



Trh 



: CA^Tc, 



(13) 



where 



C 



_3_y/^ 0.138 
0.006. 



In A 



0.01 



h/rj 



0.10 



1.5 



1.0 



0.5 



0.0 













1 — - 





0.01 



0.10 

cr/'^crl 



1.00 



Figure 2. The dependence of the dimensionless escape rate t; and expansion 
rate x on Tcr/Tcri (equations i 1 1 l and[T2l respectively). The values of x and 
5 are given relative to the efficiency of energy conduction, The top axis 
labels corresponding values of rii/rj using equation )15t and ('rh/rj)i = 
0.145. 



The approximation in the last step is appro ximately correct if we 
take A = 0.02Af dCiersz & Heggi3ll996h with = 10^ We 
turn next to the crossing time Tcr, which is defined in terms of 
the mean density inside the half-mass radius in Table [T] Since 
Tcr OC (rjJ/Af)^''^ and M oc rj, it follows that Tcr is related to 
rh/rj and its value in the evaporation-dominated regime as 



Tcrl 



[rh/r 



J 1 



3/2 



(15) 



At this point it is also convenient to give an explicit formula for 
Tcrl, which is 

3/2 



Tcrl ~ 



3GM 



(16) 



(14) 



Table[T]also introduces the 'crossing time at the Jacobi radius' , 
which may be expressed as 
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N/N„ 

Figure 3. Approximate relation between Tcr and A'^, botli normalised to 
their maximum values (equation ll9t . Evolution is from right to left. 



47rrj 



Because we now have constrained the functional forms of x 
and ^ and we have defined the time-scales involved we can solve 
for the evolution of Tci{N) and the time-dependence of both TV and 
Tcr, and we proceed to do this in the next subsection. 



2.4 Motion in the — N plane 

The evolution of N and Tcr with respect to one another can be 
found by combining the descriptions of A*' (equation I lit and fcr 
(equationll2b as derived in Section [Z2l This gives us the amazingly 
simple differential equation 

After separating the variables and integrating we find 



Tcrl 



' N' 


5/2\ 







(19) 



where iVo is a constant. Strictly, it is the value of N corresponding 
to Tcr = 0. Our model, however, is intended to describe the phase of 
balanced evolution, which begins when the cluster is compact, and 
Tcr is small but non-zero. But it is easy to see from equation il9\ 
that, when Tcr is small, A'o — A'^ ~ (2/5)(rcr/Tcri)A*'o. Therefore 
the value of A'^ at the start of balanced evolution differs little from 
A'o provided that the value of Ta then is much smaller than Tcri. 
Fig.[3]displays the relation between Tcr and A'^. The arrows indicate 
the direction of the motion along the track. The arrows are placed 
at constant intervals of time and this time dependence will be dis- 
cussed in Section|23] 



2.5 Time dependence 

Here we complete the solution of equations ([8} and The most 
compact expression of our results will involve the total evolution in 
time, from the moment when Tcr — until the dissolution of the 



cluster when TV = 0. Realistic applications of our model, however, 
require further consideration of this time scale, and so we take the 
further development of the model in two stages. 



2.5.1 In units of the evaporation time 

The time dependence of A'^ and Tcr can be solved by substituting 
the expression for r^h (equation [T3]l in the expression for A'^ (equa- 
tion [TTll, such that we find 



N 



(3/5)C 



(20) 



Crcrl 

which is constant for a given galactic environment. Thus we get the 
simple and familiar expression for the time evolution of A'^ 



TV(t) 



A^o + Nt, 
No (l 



t 

TcvO 



(21) 
(22) 



where Tovo is the total life time, which is the time it takes to reduce 
the number of stars from TVq to zero. This can be expressed as 



TevO 



-No/N, 

Trh 



(3/5)C' 



(23) 



where frh is the relaxation time of a tidally-limited cluster with TVq 
stars; i.e. we substitute N — No and Tcr = Tcri, giving 



Trh = CTVoTcrl. 



(24) 



The time evolution of Tcr is now retrieved from the expression 
for N(t) given above and the expression for Tcr (TV) (equation|19l>, 
i.e. 



Tcr(t) 



Tcrl 





t 


5/2\ 


1 - 








T cvO 





(25) 



A more familiar quantity, however, is rh, whose evolution in time 
simply follows from the combination of equations dlSI l, i22) and 
(US, i.e. 



( TcvO ) \ 


1 - 


t 






TcvO. 



5/2 



2/3 



(26) 



with fh = rjo{rh/rj)i. Here rjo is the Jacobi radius of a cluster 
with TVo stars, and rh can be thought of as the half-mass radius 
the cluster would have had if it still contained TVq stars and had 

Tcr — Tcrl ■ 

Finally, the time evolution of Trh (equation [T3]l is 



Trh(i) = Trh ( 1 — 

T cvO 



1 



t 

TcvO 



-1 5/2 



(27) 



with Trh given in equation l l24t . 

The behaviour of Tci{t), N{t),ri^{t) and Trh(t) is shown in 
Fig. in The point where expansion and evaporation are equally im- 
portant, i.e. 4 = X, happens when t ~ 0.4tcvo- (This is also the 
time at which the maximum value of Trh is reached, as Trh oc Ntct, 
if we neglect the variation of the Coulomb logarithm.) Therefore 
clusters are in an expansion dominated phase for the first 40% of 
their life, and in an evaporation dominated phase in the following 
60%. In Section ID we will estimate the fraction of globular clus- 
ters in the Milky Way that is in the expansion dominated phase. 
Similarly, we can find the turning point in the evolution of Th from 
solving rh = which gives us i(rh = 0) ~ 0.5tcvo- Therefore rh 
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shall estimate by comparison with numerical experiments of vari- 
ous kinds. 

In isolated models we neglect escape, and then it is easy to 
sh ow from equation s of Section [211 that (rh/rh)rrh — (2/3)x — 
^. I Kim et al] ( Il998h estimated the value of this quantity for iso- 
lated two-component clusters with stars of masses mi and m2, and 
found, in our notati on, that ( is higher for larger ratios m2/m\. 
More quantitatively. iGieles et alJdiOlQl use direct TV-body integra- 
tions to study the expansion of isolated clusters with various mass 
functions and they find that the dimensionless rate of expansion is 
approximately proportional to /i^'^^, where /i = mmax/jTiniin. For 
^ = 10, appropriate for globular clusters, they find x ~ 0.3, i.e. 
C ~ 0.2. 

The result that the efficiency of heat conduction (i.e. Q is 
larger in multi-mass systems is confirmed by studies on the es- 
cape rate of tidally-limited clusters. Here we neglect ^"d the 
equations of Section |2T| imply t hat the dimensionless es cape-rate 
is (Af/7V)r,.h ~ ~ -(3/5l(:. lLe"e & Goodmanl ( ll995l) consider 
various mass functions in their Fokker-Planck simulations of dis- 
solving clusters and study the effect on the dimensionless escape 
rate. Strictly, however, their results relate to the rate of escape of 
mass, i.e. to [M / M)T-rh, but we neglect the distinction here. Then 
for fJ. — 7 they find that ^, and thus also ^, is a factor of ~2 higher 
than for the equal-mass case. If we assume that ^ oc then for 
/i = 10, we again find ^ ~ 0.2, a value which we adopt hence- 
forth. It is higher than for Henon's equal-mass models by a factor 
of ~ 2.5. 

Rather than estimating (rh/rj)i directly, we consider the ratio 
C/(ni/»".i)i''^. This can be done by computing the quantity Nt^^ 
for it follows from equations l |16l l and l |20l l that 



Figure 4. Evolution of the parameters , Tct, r^^ and TjIi of clusters in a 
tidal field. Time is normalised to the total evaporation time t^vo and the 
normalisation constants on the i/-scale are discussed in the text. The dashed 
lines in the bottom two panels show the extrapolation of the early evolution 
assuming that N remains constant. The dotted lines show how rj^ and T^h 
would have evolved if the evolution would have been homologous through- 
out: rh/rj = (rh/rj)! = constant. 



increases in the first half of the cluster's life and decreases in the 
second half. 

In the expansion phase, i.e. for t <^ t^vo, it follows from equa- 
tion {27} that 



Trh(i) 



TcvO 



t 



(28) 



where we have used equations J23b and J24b . Hence we arrive at the 
interesting conclusion that the relaxation time of all such clusters 
is the same at a given time. Furthermore, since this relationship is 
represented by the dashed line in the lowest panel of Fig. [4] it is 
clear that this gives an upper limit on the relaxation time at time t 
for all clusters satisfying the assumptions of our model. We return 
to this point in Section [33] 



2.5.2 In units of physical time 

To describe the evolution in units of physical time we need to know 
the value of Tcvo, which by equations l ll6t and l l23t depends essen- 
tially on C,, (rh/rj)i and the ratio rj/M (in addition to A'^o). The 
last of these three quantities is set by the tidal field of the Galaxy, 
while the first two parameters will be taken as constants, which we 



(3/5)C 



V2C(rh/r-,T)f ' ■ 
71 C 



(rh/rj 



(29) 
(30) 



where we used C — 0.006 in the last step (eauationlI4l(. This is the 
escape rate corrected for the tidal density, and should thus be the 
same for clusters at different locations in a galaxy. In Table [2] we 
summarise a number of results on this quantity from the literature, 
expressed in this way. 

For Henon's equal-mass model ~ 0.0743 and (rh/rj)i = 

-95. The Fokker- 



0.145, Sections O and [Ta we find Nt^^ 

Planck models of equal-mass clusters of iLee & OstrikeJ h987l) 
start relatively compact, and so nearly the entire evolution of their 
clusters is spent in the post-collapse phase. Their results imply that 
Nt^,. — —76. On the other hand their values for ^ or (i.e. the es- 
cape rate scaled by the half-mass relaxation time, as in equation {8]( 
are slightly larger than Henon's value. But from their Fig. 2 we see 
that in their models (rh/rj)i ~ 0.2, which is slightly larger than 
the value (rh/r,j)i — 0.145 in the model of Henon. This explains 
their slightly lower value Nr^r (equationl30t. It also illustrates that 
5 is not the only parameter that determines the total life time in 
physical units. 

Note that the results from studies based on A'^-body simula- 
tions in Table|2]have an additional factor proportional to N^^'^. This 
A'^-dependence of the escape rate is reasona bly well understo od, in 
terms of the time taken for stars to escape l lBaumgardtlboOlh . and 
is included in the more refined models in Appendices IaI and ICl 
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Table 2. Overview of literature results for the escape rate Nt^^ in post- 
collapse/balanced evolution for clusters with different ^ = rrimax / mm- 
Values are de rived from the results of (l llHenon|jT96l|) . (2) Lee & Ost rikeJ 
<1987l). (3) iLee. Fahlman & RicheJ I199ll). (4) iHeggie et al UgsiT 
rSl lGieles & Baumgardj <2008ll . f6l lLamers. Baumgardt & Gielej fcOld) . 



Reference 


* ' cr 






Comments 


(1) 


95 




1 


Theory 


(2) 


76 




1 


Fokker-Planck 


(3) 


200-250 




7 


Fokker-Planck 


(4) 


290 {N/IQ^ 


jl/4 


15 


Af-body 


(5) 


300 (TV/IO^ 


1/4 


30 


Af-body 


(6) 


270 (iV/10^ 


1/4 


~10 


A'^-body+stellar evolution 



Notes on the data. 



(2) Based on their value of A (p. 128). 

(3) Based on their value of rov (p. 458) and a mass-function index x = 1. 

(4) Based on a run with N = 65 536. The dependence on A'^ has been 
assumed from lBaumg^ bOOlh . 

(5) Based on the results of their Fig. 2. This study found a slightly 
stronger A'-dependence; the N^^^ was assumed and matched to the runs 
with the largest number of particles (N = 32 768). 

(6) From a fit to the data presented in their Fig. 4 (Lamers 2010, private 
communication). This study also finds a slightly stronger Af -dependence; 
the Af ^/^ was again assumed and matched to the runs with the largest num- 
ber of particles (A^ = 131 072). 
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Figure 5. Isochrones (full lines) and evolutionary tracks (dashed lines) ex- 
tracted from the model described in Section [2] applied to clusters orbiting 
at a radius of /Jq = 8 kpc in an isothermal halo with circular velocity of 
Vc = 220 km s~^. Arrows along the tracks indicate the direction of evolu- 
tion and have a length of 0.2 dex in time. The central isochrones has an age 
of a Hubble time in = 13 Gyr. 



3 COMPARISON TO THE GLOBULAR CLUSTERS IN 
THE MILKY WAY 

3.1 Model parameters 

In this section we compare the basic model of the previous sec- 
tion to parameters of Milky Way globular clusters. We convert 
the model results for A'^ to M assuming m = 0.5 and convert 
Tcr to ph (equation [6}. The tidal field parameters Tcri and tI^ 
(equations 1161 and 117b are converted to Rg by approximating the 
Milky Way potential by that of an isothermal sphere with a (con- 
stant) circular velocity of Vc ~ 220kms~^ and using Henon's 
value (rh/rj)i = 0.145. The resulting formulae are given in Ap- 
pendix|B] With these parameters fixed, the speed of evolution is set 
by the value of C, which is chosen to be = 0.2; this is appropriate 
for clusters with a globular cluster-type mass spectrum (/x ~ 10, 
Section l2.5.2t . As a consistency check, by equation l l3Qt , we derive 
NtI, ~ -256, which may be expressed as (see Appendix iBt 



MRg ~ -20 Mq Myr"^ kpc. 



(31) 



The first of these is in reasonable agreement with the evaporation 
rates found in both A^-body and Fokker-Planck models with a com- 
parable mass function (TableO. 

From these data we can obtain the value of Tcvo from substi- 
tution into the equations of Section 12.5.11 This then gives us the 
relation between ph, M and Rq at a given time (isochrones) and 
the evolution of the individual parameters in time ph(Afo, Rojt), 
M{Mo, RG,t) (tracks). In Appendix|B]we provide these relations 
and also some of the intermediate steps and additional relations for 
rh(A/, Rg) and rh(Mo, /?g, i) that are not shown here. 

In Fig. [5] we illustrate the behaviour of the evolutionary 
tracks (dashed lines) together with various isochrones (full lines) 
in these more appealing quantities for clusters at Rg ~ 8 kpc. 
The isochrones span 2 dex in age centred around a Hubble time 
(tu — 13 Gyr). The asymptotic behaviour of the isochrones at 



the extremes is easy to understand. In the expansion dominated 
phase (right-hand side of Fig. |5]( all clusters have evolved to the 
same relaxation time (cf. equation I28t and thus ph oc (M/t)^. 
The proportionality, i.e. the absolute vertical positioning of the 
lines, is set by the value of In the evaporation-dominated regime 
(left-hand side of Fig. [5j ph has adjusted to the tidal density and 
ph oc Rq^. This is because clusters are in the homologous phase 
where r^/r] — (rh/r,j)i and therefore ph oc p.j; while for an 
isothermal halo the lacobi density, pj scales with Rg as pj oc Rq^ 
(Appendix iBt. The tracks (dashed lines) show that clusters initially 
move down in this diagram, because they are expanding without 
losing much mass. Evaporation becomes more important than ex- 
pansion when ph approaches its minimum value and the clusters 
start moving to the left. In the limit of t — >■ cxd both the tracks 
and the isochrones are horizontal lines. The speed of the horizon- 
tal motion is set by the value of and Rg because it follows from 
equation ( |20b that M (x (/ Rg ■ 



3.2 The evolutionary state of Millsy Way globular clusters: 
The relative importance of expansion and evaporation 

In order to see if these types of predictions are relevant for real 
globular clusters we will compare our results to the globular clus- 
ters of the Milky Way. We use the 2003 version of the lHarrisl ( [l996l) 
catalogue. It contains entries for 150 globular clusters, and for 141 
of them a luminosity, radius and galacto-centric radius determina- 
tion are available. To convert luminosity to mass we a dopt a mass- 
to-light ratio of 2 jMcLaughlin & van der Marell2o"ol) and we mul- 
tiply the projected half-light radius by 4/3 to correct for the effect 
of projection (Spitzer 1987) and get an estimate for r^. Note that 
by doing this we assume that light traces mass, which is not nec- 



essarily true if the cluster is ma ss segregated (e.g. lLee et al.lll99ll ; 
lHurlevil2007h . iBonatto & Bic j (2008) study this effect in a sam- 
ple of 1 1 globular clusters. They compare the half-light radius to 
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Figure 6. Isochrones based on the balanced evolution models of Section [2] In the left p anel 5 isochro nes for different Rq values are shown and in the right 
panel 5 isochrones for different M. Both panels show the 141 globular clusters in the iHarri^ (l9%) catalogue from which M, and Rq estimates can 
be derived. The filled circles are the 48 clusters classified as being in the evaporation-dominated regime according to equation (32). The traingles are the 
remaining 93 classified as being in the expansion dominated regime. 



the radius containing half the n umber of stars and find that the lat- 
ter is typically ~ 10% larger. iHurlevI j2007h finds from A'^-body 
simulations that rh can be about a factor two larger than the pro- 
jected half-light radius, due to mass segregation. This implies that 
for some globular clusters rh could be a factor 2/(4/3) = 1.5 
higher than what we derived in Section [T2l This might affect clus- 
ters of different masses in different ways; this is further discussed in 
Section |43] For now we assume that light traces mass and that the 
correction for projection provides a sufficiently accurate estimate 
of the 3-dimensional half-mass radius rh. 

The first thing we determine from the data is the fraction of 
globular clusters that are in the expansion dominated phase. We 
adopt the definition of the boundary between these two phases as 
the moment where the time derivative of Trh is zero as outlined 
in Section [2. 5. II In terms of Tcvo this point is when the cluster's 
age is 0.4revo approximately. This means that we need to find 
the parameters of clusters that at present have a life-expectancy of 
Tev — l-5iH, i.e. with 60% of their evolution still ahead in time. 
Note that our estimates are based on an assumption that all clusters 
have spent the last at the Rq where they are now and that the 
potential of the Milky Way halo has been constant. This is not gen- 
erally expected to be true and the consequences of this assumption 
are discussed in Sections 14.21 and 14.31 The exercise is, therefore, 
merely intended as a first order approximation to see if it is at all 
reasonable to expect that globular clusters are still in the expansion 
phase. 

With the adopted model parameters (Section l3.lt and tn = 
13 Gyr we find that clusters with 

M>1O^M0(^^^ (32) 

should still be in the expansion-dominated phase. (Equation 13 II is 
the easiest route to this result.) This relation is satisfied by 93 of the 
141 clusters. It follows that the remaining 48 clusters are in the 



evaporation-dominated phase. This perhaps surprising result has 
some interesting consequences. The most important one is that the 
present day densities of the majority of the globular clusters fol- 
low (roughly) from the self-similar expansion model for isolated 
clusters: ph oc . Here we have assumed that all clusters have 
the same age. This scaling relation should be universal, because it 
is driven by internal 2-body relaxation; therefore a similar scaling, 
with the same proportionality, should also hold for extra-galactic 
clusters. Moreover, in extra-galactic cluster samples the fraction of 
clusters in the expansion-dominated phase is probably larger; they 
are easier to detect because they have (on average) higher mass 
(equation 1 3 2b. 

1/2 

The prediction that a oc M scaling must hold for the ma- 
jority of the Milky Way globular clusters is one of the main results 
of this paper. More detailed discussion on its implications and con- 
sequences are postponed to Section|4] First we show in Section [33] 
that this scaling is indeed found for globular clusters and that the 
proportionality agrees with the parameters of energy conduction 
derived in Section|2] 

3.3 Isochrones and Milky Way globular clusters 

Because all globular clusters have roughly the same age we focus 
on isochrones with an age of — 13 Gyr, rather than the evolu- 
tionary tracks. In Fig.|6]we show the isochrones in ph vs.M (left) 
and ph vs.Rg (right) diagrams together with the 141 globular clus- 
ters for which data are available in the Harris catalogue. The clus- 
ters that are in the evaporation-dominated phase are shown as filled 
circles and the clusters that are still expanding are shown as trian- 
gles. 

In the left panel isochrones for clusters at different Rq be- 
tween 1 kpc and 100 kpc are shown. The isochrones roughly en- 
compass the data. The 100 kpc isochrone clearly shows the asymp- 
totic pI/'^ oc M behaviour following from expansion, which 
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Figure 7. Similar to left panel of Fig.|5] but now the data is sorted in Rq 
and sub-divided in four (roughly) equal samples. The range in Rq values 
is indicated in each panel. The full lines show isochrone predictions for the 
median Rq values in the samples and the dashed lines show the isochrones 
for the minimum and maximum Rq values. 



Figure 8. Similar to right panel of Fig. [5] but now the data is sorted in 
M and sub-divided in four (roughly) equal samples. The range in AI val- 
ues is indicated in each panel. The full lines show isochrones predictions 
for the median values of AI in the samples and the dashed lines show the 
isochrones for the minimum and maximum values. 



roughly follows the lower envelope of data points. In the outer 
halo the tidal field is so weak that all clusters with M > 10'* M© 
have not expanded up to their tidal boundary yet. In the right panel 
the densities are shown as a function of Rq together with five 
isochrones for different masses. As was the case for the isochrones 
in the ph vs. M plot, these isochrones also roughly encompass the 
data. The asymptotic behaviour of the isochrones in both diagrams 
is given by labels in the two diagrams. 

To further illustrate the comparison between data and theory 
we subdivided the Milky Way sample in four samples containing 
roughly equal numbers, but sorted in M and Rq. We then compute 
the isochrones based on the median, the minimum and maximum 
value in each sample. The results for different Rq bins and different 
M bins are shown in Fig. |7] and Fig. [8] respectively. The general 
relation between M, pi^ and Rq is nicely matched by the model. 

Our model is based on the assumption of balanced evolu- 
tion (Section [TJ, which led to the conclusion that ;$ 0.3tH 
(Section 12.5.11 equation I28t . A logical conclusion of our results, 
then, is that the evolution of objects with r^h ^ 0.3tH is unbal- 
anced, presumably implying that they have formed with a long 
Trh (this scenario was e xplored using nume rical simulations by 
iHurlev & Mackevll2010l : [Zonoozi et alJIloT ll). If we draw a hard 
line at this limit then we find 23 globular clusters for which this 
holds. However, many of these clusters are fairly close to the 
Trh = O.StH line (see left panel of Fig. |6j. The majority of the 
points which are substantially below the lowest isochrone have a 
low mass, and given the uncertainties in radius, mass and distance 
estimates for these objects we should not exclude the possibility 
that their positioning in the diagram is consistent with the model 
within measurement uncertainties. If these objects are in balanced 
evolution, they have undergone a lot of dynamical evolution, and 
could therefore be mass segregated, and perhaps estimates for A*', 
m and rh are systematically affected by this. In Section |4] we fur- 



ther discuss this issue of the clusters with Trh > O.Sfn- Obvious ex- 
ceptions (i.e. high-mass clusters well below the lowest isochrone) 
are NGC 5139 (cjCen) and NGC 2419 with ~ 16 Gyr and 
Trh — 54 Gyr, respectively. These are the only two objects in the 
Harris catalogues with Trh > tn- Due to their high mass it is un- 
likely that the data of these objects suffer from measurement un- 
certainties and it is more likely that these objects have properties 
resembling their initial conditions (besides some adiabatic e xpan- 
sion because of mass-loss due to stellar evolution perhaps, iHillsl 
.198a : .Gieles etal..201Q) . 



4 SUMMARY AND DISCUSSION 

In this study we provide a description of the evolution of mass 
and half-mass density (or radius) of (globular) clu sters. The two 
self-similar solutions provided bv .Henon, 4l965h and lHenol] ( Il96lh 
for clusters evolving in isolation and in a tidal field, respec- 
tively, are combined into a pair of Unified Equations of Evolution 
(UEE), which yield analytically continuous evolutionary tracks and 
isochrones. The key ingredient is that clusters are assumed to con- 
duct a constant fraction of their total energy per half-mass relax- 
ation time-scale. This is assumption is justified by the fact that both 
Henon models, i.e. those that describe the extremes ('Henon"l96l|, 
1965), have a very similar {E/\E\)T^t, ~ 0.08 (see Section ITTTl. 
The central energy source could be assumed to be hard binaries in 
the centre of the cluster, but the details of the energy production are 
not of interest to us, since we assume that the evolution is balanced, 
i.e. that the amount of energy produced by the core is determined by 
how much energy can be conducted through the half-mass bound- 
ary. This important view on the evolution of c ollisional gravita- 
tional systems was introduced bv lHenonI j 1 9751) and can be com- 
pared to Eddington's view on stellar energy sources. The bottleneck 
in the energy emission from a star is the radiative transfer through 
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the envelope. The luminosity of the star is thus set by how much ra- 
diation can be transported through the envelope and the nuclear fu- 
sio n rate in the core adjusts to this (see lLvnden- Bell &Woodl 19681 
and llnag^i & Lvnde n-Bell 1981 for more details on this analogy). 

The efficiency of energy transport in a cluster depends on the 
mass spectrum of the stars. For a mass function typical of a glob- 
ular cluster it is about a factor of ~ 2.5 higher than for the equal 
mass clusters considered by Henon. It is this efficiency that deter- 
mines the positioning of the isochrones and the speed with which 
clusters move along the tracks. The solutions are compared to the 
parameters of Milky Way globular clusters and reasonable agree- 
ment is found for almost all well-observed clusters. We interpret 
this agreement as evidence for the fact that the evolution of almost 
all globular clusters is balanced, meaning that there is a central en- 
ergy source that supplies the energy needed by 2-body relaxation. 
Exceptions are NGC 2419 and cjCen because their relaxation time- 
scales are much longer than a Hubble time. The evolution of these 
clusters is unbalanced. We find that roughly two-thirds of the Milky 
Way globular clusters are in the expansion dominated phase, i.e. 
their evolution is not yet seriously affected by the Galactic tidal 
field. These clusters align along lines of constant relaxation time 
implying that ph oc KP . This idea is further supported by the prop- 
erties of intra-cluster globular clusters in the Virgo Cluster, which 
have no tidal boundar y. These have masses of M ~ lO'^ Mq and 
radii of a few parsecs jWilliams et"^l2007h , comparable to Milky 
Way globular clusters. 



4.1 Are all globular clusters in post-collapse evolution? 

We have referred to balanced evolution when we talked about clus- 
ters with a central energy source. For clusters without primordial 
binaries and stellar evolution this is generally referred to as post- 
collapse evolution, the energy being provided by binaries formed 
in three-body interactions. It is interesting to enquire whether these 
concepts are also related in more realistic models. 

At first sight the proportion of Galactic globular clusters which 
are classified as post-collapse objects is quite inconsistent with our 
conclusion that the evolution of almost all globular clusters is bal- 
anced. Only ~ 20% of the globul ar clusters possess the stee p cusp 
in their surface-brightness profile jPiorgovski & King|l986h that is 
predicted by models of post-collapse evolution. We should bear in 
mind, however, that such predictions are based on studies of equal- 
mass systems. In a system with an evolved stellar population, how- 
ever, it is possible that the distribution of the remnant population is 
steeply cusped, and that the central distribution of those stars which 
dominate the surface-brightness exhibits only a shallow cusp. This 
has been shown in many s tu dies since a t least the work on M15 by 
llllingworth & Kind ( ll977h ; lGoodmanl h984h gives a general and 
simple theoretical explanation. 

Several independent arguments support the view that the 
surface-brightness profile is not a reliable guide to the post- or pre- 
collapse status of a cluster. Cohn & Hut ( 1984) argued that globular 
clusters with central relaxation times less than 10* yr have already 
undergone core-collapse, which is clo se to half of thei r samp le of 
146 Galactic globular clusters. Next, IPe Marchi et alj ( |2007|) . on 
the basis of a surprising correlation between the concentration of a 
cluster and the slope of the main sequence mass function, argued 
that the number of post-collapse globular clusters has been seri- 
ously underestimated. Finally, specific examples are being revealed 
by detailed modelling of individual clusters. A good example is the 
cluster M4, whose surface-brightness profile is well described by a 



classical King model, but dynami cal modelling suggests that it has 
already undergone core-collapse ( iHeggie & Giers 320081) . 

While these arguments reduce the gap between the propor- 
tion of post-collapse clusters, on the one hand, and the propor- 
tion exhibiting balanced evolution on the other, they do not imply 
that the two concepts a re identical. In another modelling exercise 
lGiersz&Heggia ( l20Tlh argue that 47 Tuc is far from core collapse, 
and yet the model exhibits what we would call balanced evolution, 
the required energy being supplied by mass loss from stellar evolu- 
tion and, to a lesser extent, heating by primordial binaries. In this 
model, what is called 'core collapse' results from the long-term ex- 
haustion of these energy sources. 

Models of globular clusters show an early rapid collapse of the 
core radius (e.g. Fig[Tll caused by mass segregation. It is also exhib- 
ited by the model of 47 Tuc which we have just been discussing. 
This phenomenon is the closest analogue of the classical core col- 
lapse exhibited by single-component models, and perhaps the main 
lesson of this discussion is that we need a better, agreed definition 
of what we mean by core collapse. 

4.2 Cluster orbits 

In our model we have made the implicit assumption that cluster 
orbits ar e circular in the M ilky Way potential. This is of course 
not true. IPinescu et al. I ( fl999) determined the orbits of 38 globu- 
lar clusters and they find a median eccentricity of e ~ 0.65. The 
question is to what extent this affects our interpretation. 

Clusters that are in the expansion-dominated phase are not af- 
fected by the tide, and because we have shown that this is the case 
for the majority of the globular clusters in the Milky Way, the de- 
tails of the shape of the orbits should not affect the result that the 
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empirically established scaling of oc A4" is a result of expan- 
sion. However, our estimate of the fraction of clusters in the ex- 
pansion dominated phase was based on the assumption of circular 
orbits. This fraction could very well be smaller because a cluster 
that is in the expansion phase based on its density and its current 
galacto-centric position could be in the mass loss dominated phase 
at peri-centre (-Rp) when we allow for the possibility that e > 0. 

Despite what is commonly assumed, conditions at Rp are not 
relevant, either for the struc ture of the clu ster or for the evapora- 
tion rate. On the first point. iKtipper et al] (2010) showed that the 
density profile of the bulk of the stars adjusts to the mean tidal 
radius along the orbit, rather than the tidal radius at perigalacti- 
con. On the second point, again through use of A'^-body simula- 
tions, it has also been shown that the escape rate of a cluster on an 
oval galactic orbit is close to that of a cluster on a circular orbit 
with a radius intermediate betwee n Rp and the apogalactic dis- 
tance l lBaumgardt & Makindl2003l) . Though it has not been shown 
that the galactocentric radius determining the evaporation coincides 
with that characterising the structure, the systematic error in fitting 
the evolution of a cluster to its present radius must be considerably 
smaller than would be the case if structure and mass loss were set 
at Rp. Nevertheless, whatever the magnitude of this effect, it does 
not help to explain those clusters which, in Fig|6](left panel), lie 
below the envelope of the isochrones. 

4.3 Static galaxy potential 

We have assumed that clusters have been evolving for a Hubble 
time at the galactocentric position where they are observed now. 
This is an oversimplification, because in the hierarchical merging 
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scenario clusters can lorm in dwart galaxies that are later accreted 
in larger haloes jKravtsov & Gnedi J2OO5I : IPrieto & Gnedinll2008l : 
iMackev et al.ll201(]h . Most of the galaxies are in place at a redshift 
2: ~ 2 — 3, hence the assumption of a static galactic potential only 
affects the clusters that were brought into haloes in late accretion 
events. In fact, there is the potential of using the densities of clusters 
to say something about the accretion history. Imagine a cluster in 
a dwarf galaxy with a strong tidal field. Suppose that in this host it 
has already reached the evaporation-dominated phase, i.e. its den- 
sity has adjusted to the tidal density in the dwarf galaxy (which 
is the case for the clusters as sociated with the Sagittarius dwarf 
galaxv, |Penarrubia et alj|200^ . If the dwarf galaxy is accreted in a 
larger host, such as the Milky Way, the new tidal field that the clus- 
ter experiences is in most cases weaker than in the previous host. 
The cluster will start expanding again, but its density is higher than 
the prediction based on a Hubble time of evolution at that Rg . A 
target to which this might apply is Palomar 1 (see the discussion in 
iNiederste-Ostholt et al.ll20IOi) . 



4.4 Model validity at younger ages 

Throughout this study we have assumed a constant value of ~ 0.2 
for all ages. Gieles et al. ( 2010) found values for C, as large as 2 for 
clusters with a full mass spectrum (/i = 1000), which would be ap- 
propriate for very young star clusters (< 3Myr). We have also as- 
sumed that the evolution is balanced at all ages. Therefore we need 
to bear in mind that the results presented here should not be inter- 
preted in a literal way for clusters with ages much younger than a 
Hubble time. The results presented here underestimate the speed of 
balanced evolution for younger clusters. For clusters that have not 
yet reached the balanced evolution phase our results overestimate 
the speed of evolution. The expansion for clusters with young ages 
is considered in more detail by Giele s et al.. ( 2010) and we refer to 
that paper for an analytic description of the radius as a function of 
initial r^h considering the effect of the changing stellar mass func- 
tion and the transition from un-balanced to balanced evolution. A 
time-dependent functional form for C, co uld in princip l e be f ound 
from the time-dependent parameter Xt of iGieles et al.l ( l2010h . We 
refrain from doing this, however, because Xt depends on the de- 
tails of the stellar mass function, which is affected by both stellar 
evolution and by the loss of low-mass stars over the tidal boundary. 
The l atter was not taken into account in the models of jGieles et al.l 
I2OIOI) . We therefore leave these details of the full evolution for a 
future study. 



4.5 The stellar mass function 

In Section [331 we found that there are a number of low-mass clus- 
ters with O.Stn '''rh J; ^h, i.e. their relaxation times seem in- 
consistent with the predictions for expansion with C, — 0.2. Several 
Palomar clusters are in this group. We now consider a number of 
factors which might account for such cases. 

Because these clusters are faint, they could be close to final 
dissolution and m ight have stellar mass functio ns that are depleted 
in low mass stars ( iBaumgardt & Makindl2003h . We assumed (Sec- 
tion |2.5.2"l ( that the value of which controls the rate of evolution, 
depends only on the range of the mass function, i.e. the ratio /i of 
the maximum to minimu m mass, and this doe s not depend on the 
slope of the mass function. lGurkan et al.l ( l2004h showed that for the 
rate of core collapse it is not n = mmax/'Timin that sets the speed 
of evolution, but rather mmax/'Ti, where m is the mean mass. If 



the same consideration applies to balanced evolution, then for flat- 
ter mass spectra, i.e. lower mma^/fh, we would expect a slower dy- 
namical evolution, i.e. smaller C,. By equation ( |23t this increases the 
ratio Tcvo /Tcri. By use of equation ( I24t we see from equation ( 128b 
that the effect is to reduce the ratio Trh/t in the expansion phase. 
Thus for a given mass, clusters should be denser at a given age. On 
the other hand Fig|6](left panel) shows that the clusters of concern 
are less dense than expected. We conclude, therefore, that a flatter 
mass function can not be responsible for an additional decrease of 
the densities of these low mass clusters. 

Another consideration is the A'^-dependence in the escape rate 
which is implicit in the A'^-body results quoted in Table [2] This 
affects the expansion rate x of clusters that are well into the 
evaporation-dominated regime (Appendix [Aj and therefore does 
not help us to understand the discrepant objects. 

We propose two alternative speculative explanations. Firstly, 
it could be that these clusters had a high retention fraction of black 
holes and neutron stars. If the most massive stars are black holes 
with m ~ 10 — 20 M©, the ratio ^ would be higher and the dy- 
namical evo l ution would be faster. As mentioned in Section |2.5.2| 
iGieles et al.l ( I2OIO I) show that C, oc /x^^^ approximately. For a mass 
function between 0.1 Mq and 1 M0 we have adopted C, — 0.2. For 
a cluster with a few additional black holes with masses of 10 Mq 
we would have (" approximately a factor of ~ 3 higher. This would 
help to make all of these clusters consistent with balanced evolu- 
tion. 

Secondly, we have assumed that the mass-to-light ratio is the 
same for all clusters. But if the cluster is depleted in low-mass stars 
we surely overestimate and underestimate m by using a canon- 
ical mass-to-light ratio. Because r^h oc (N/rh)^^'^ we thus over- 
estimate Trh. Thus if, according to our estimates, a cluster has a 
half-mass relaxation time-scale > O.Stn, it may be that its actual 
Trh is more consistent with our theory. It would be interesting to 
consider these effects in more detail with numerical simulations. 

Lastly, estimates of Trh are very sensitive to uncertainties 
in the distance to the cluster, D. The radius in pc depends lin- 
early on D and the mass in solar masses quadratically, such that 
Trh oc D^^^. Because ph oc D"^ and M oc the scatter of data 
points due to uncertainties in D is orthogonal to lines of constant 
relaxation time (left panel of Fig.[6]and Fig.|7]l. 



4.6 Comparison to other work 



The relation of globular clusters to lines of constant relax- 

1/2 

ation time (pj^ oc M) has been noted for Milky Way 
globular clusters and for extra-gala ctic globular cluste rs (e.g. 
Gnedin & Ostri ke3ll997l: iBarmbv et a l. 2007; McLaughlin & Falll 
20081 : iGeorgiev et al.l |2009|) . It is usually given a different in- 



terpretation, however. Because for homologous clusters the life- 
expectancy is a constant number of relaxation times (1/0 the lines 
of constant Trh have been interpreted as boundaries to the expected 
distribution of surviving clusters. The first to provide this interpre- 
tation were probably Fall & ReesI d 19771) who showed the depen- 
dence on mass and radius of various disruptive agents. Given that 
the time-scales for relaxation driven evaporation (of homologous 
clusters), disruption by tidal shocks and dynamical friction all de- 
pend in a different way on cluster mass and radius one can con- 
struct a 'survival' or 'vi tal' triangle in an r^^vs.M diagram (e.g. 
iGnedin & OstrikeJl997l) . or alternatively phVS.M. Clusters within 
this triangle are likely to survive for another Hubble time. 

In this paper we have given a different interpretation to the re- 
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lation between cluster data and lines of constant Trh- Because the 
majority of the clusters have not reached the homologous (tidally 
limited) evaporation-dominated phase yet, the alignment of the 
globular cluster data for these clusters with the — O.Stn rela- 
tion results from expansion driven by 2-body relaxation when evo- 
lution is balanced. There would be no globular clusters with shorter 
relaxation times if there was no tidal field that can stop the ex- 
pansion. In reality clusters can not expand indefinitely, and the re- 
laxation time starts decreasing once evaporation starts to dominate 
over expansion. This is why there are clusters with Trh < O.Stn. 
Note that this is just the opposite of the usual interpretation, which 
implies that only clusters with longer relaxation times should be 
observed! 

Because Trh decreases linearly in time in the homologous 
phase, there should be a uniform distr ibution of Trh values betw een 
and ~ O.Stn. This was proposed bv lLee & Goodm^il j 19951) and 
they show that this indeed roughly holds. Since in the homologous 
phase Trh oc M, there should al so be a uniform distribution of M a t 
low mass, or AN /AM ~ const ( iHenonll 196 ll: iFall & Zhandl200lh . 
whi ch is indeed fo u nd in many galaxies ( Jordan et al.ll2007t) . 

'Kii pper et alj j2008l) recently presented results of A'^-body 
simulations of initial compact clusters in a tidal field. They refer 
to the phase where the half-mass radius has adjusted to the Jacobi 
radius as 'main sequence' evolution. Using again the analogy with 
stars we could, however, refer to the entire balanced evolution, i.e. 
expansion and evaporation, as main sequence evolution. The initial 
unbalanced evolution without a central energy source would then 
be the 'pre-main sequence evolution'. On t he other hand perhaps 
the an alogy should not be pushed too far; for lLvnden-Bell & Wood 
( Il968l) the end of core collapse signalled the start of the red giant 
phase] 

IBaumgardt et al.1 Jioioh identified two populations of clusters 
based on the ratio of the half-mass radius over the Jacobi radius 
{Th/rj). About half of the clusters beyond > 8kpc is strongly 
under-filling their tidal radius (i.e. rh << rj), while the other half 
have 0.1 < rh/r,j < 0.3. They find a strong correlation between 
mass and the membership of these groups: the under-filling group 
contains significantly more-massive clusters than the 'tidally fill- 
ing' group. This separation is because low-mass clusters have al- 
ready expanded towards their tidal boundary and all have a similar 
ratio rh/rj . Their more massive counterparts are still expanding to- 
wards their tidal boundary (see right panel of Fig.|6](. The densities 
of clusters in the expansion-dominated regime within Rg < 8 kpc 
are more similar to those of the clusters that are already in the 
evaporation-dominated regime (Fig. [6]( and it is, therefore, harder 
to make a distinction between under-filling and filling based on the 
empirically determined ratio rh/r,j. 



4.7 Implications for the use of globular clusters as standard 
rulers for distance estimation 

Globular cluster properties have been used as standard candles for 
distance estimates. Traditionally the peak of the globular cluster lu- 
minosity funct ion (GCLF ), because of its near universality, is used 
for this (e.g. Harrislll99ll) . However, even ;/the initial mass func- 
tion of globular clusters is universal, the detailed shape of the GCLF 
will be environment dependent after a Hubb le time of dynami- 
cal evolution (Vesp erini & Heggie 1997 ; Ostrik er & Gnedinll 19971 ; 
IVillegas et all 20 lol) . I Jordan et all l2005b suggested that the radius 
distribution can be used as a standard ruler. Alternatively, we sug- 
gest to use a combination of the two as a distance ruler; namely 
the relation between apparent luminosity and radius. This can be 



done because for the majority of clusters the mass-radius relation 
is independent of environment and second order effects such as the 
stellar mass function can in principle be taken into account (see also 
[jordaa.2004.) . 



4.8 Relation to fundamental plane relations 

With all the theory in place we are able to test if empirically 
established 'fundamental-plane relations' (e.g. Diorgovskii ri995l ; 
Burstein et al. 1997; McLaughlin 2000) can be attributed to bal- 
anced evolution. Because we do not make predictions for the cen- 
tral velocity dispersion, we will focus here on the mass-radius rela- 
tion. 

It is often stated that the ra dii of globular clus t ers d o 
not correlate with their mass (e.g. Ivan denBerghet"d]|l99lh . 
This statement is perhaps a slight oversimplification, because 
there may be a negative correlation, especially for clusters in 
the o u ter halo Ivan den Be rgh & Mackevll 20041 : IBaumgardt et al.l 
I2OIOI) . McLaughlin ( 2000) showed that, if a correction is made for 
the galacto-centric radius, the cluster half-mass radius is approxi- 
mately constant, i.e. independent of mass: rhZ-Rc* ~ constant. 

Strictly speaking the theory presented in this paper can not 
explain a radius that is perfectly independent of mass. The relation 
between mass and half-mass radius depends on the regime the clus- 
ter is in. For clusters in the expansion dominated regime we expect 
a constant relaxation time, i.e. rh oc M^^^^ independent of Ro, 
while for clusters in the mass loss regime we expect a stronger de- 
pendence on Rg than on M : rh oc M'^^^R^J^ (or rh oc M^^'^R^J^ 
when we use a slightly different escape criterion, see Appendix IaIi. 

B 

If we split the Milky Way globular cluster sample, using the 
boundary between these two extremes as described in Section [J!2l 
then we can search for these different (bivariate) relations in the 
data and see if they indeed hold. We use a straightforward linear 
regression fit to find how log rh depends on log M and log Rq ■ 
We exclude cjCen and NGC 2419 because the balanced evolution 
does not apply to these objects. We also exclude NGC 6540, with 
a density of ph — 7 x 10'* Mqpc^"', because it is an extreme out- 
lier in all diagrams. From a bivariate fit to the 47 clusters in the 
evaporation-dominated regime we find 

log(rh/pc) = -0.16+ (0.11 ±0.08) log(A4"/M0) 

+ (0.44 ± 0.11) log(iiG/kpc). (33) 

This shows a weak positive correlation with AI close to the pre- 
dicted index of 0.17(0.33) for x = 3/4(1) (see Appendices lAl and 
|C]for the significance of x) and a relatively strong dependence on 
Rg, in rather poorer agreement with the predicted index of 0.67. 
On the other hand we find for the 91 clusters in the expansion dom- 
inated regime 

log(rh/pc) = 1.82 - (0.25 ±0.05) log(M/M0) 

+(0.27 + 0.05) log(i?G/kpc). (34) 

This is consistent with a constant relaxation time (rh oc M"^'''^, 



^ The correlation between rh and M becomes weaker still if we take into 
account the Coulomb logarithm. If we approximate In A by a weak power 
of M: In A oc M°-'^^ then we find in the expansion dominated regime that 
rh oc In A^/^M"^/^ ~ M~^ '^^ and in the evaporation dominate regime 
rh oc Mi/3/lnA2/3 - M°-23 (or rh oc Mi/^/lnA^/S ~ A/"-"'' for 
X = 3/4, see Appendix lAl. 
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equation [3j, but there is still some dependence on Rg, whereas 
pure expansion would not give rise to any Rq dependence. 

The balanced evolution we consider here is thus responsible 
for a slightly complicated correlation between radius and mass. 
The theory predicts that the radius has a slight positive correlation 
with mass for M sC 10^ M0(4kpc/_RG), while for more massive 
clusters the correlation is negative and consistent with a constant 
relaxation time relation. This subtle behaviour of the mass-radius 
relation is indeed found (at least qualitatively) for the Milky Way 
globular clusters, as we have just seen. An approximately mass- 
independent radius (at a given Rg) is thus only true to first ap- 
proximation for the population as a whole. The strong correlation 
of Th with Ro for the cluster population as a whole does not seem 
consistent with our finding that the majority of the clusters are still 
expanding. However, we note that the scaling rh oc .Rq* is very 
sensitive to the handful of very low density clusters in the outer 
halo (> 50 i?G). 
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APPENDIX A: ALTERNATIVE ESCAPE CRITERION 

In Section |2] we have solved the simple case in which the escape 
rate is constant during the entire evolution. This resulted from the 
fact that we used the scaling ^ oc Tcr, cancelling out the Tcr term of 
Trh in equation jilt . The constant escape rate is to first order a good 
approximation, but a slightly improved description for ^ would be 
the exponential function mentioned in Section l2!2l In principle all 
derivations of Section|2]can be repeated using different functional 
forms for ^. Together with the constraint that C, — constant we can 
always solve for x- For the exponential function mentioned before 
we would have to proceed with numerical integrations of the vari- 
ous differential equations, which we will not do here. 

Stars that have enough energy to escape take a finite time to 
do so, and this has a profound effect on the escape rate. This time 
scale can be very long f or stars with energies only s lightly higher 
than the escape energy jFukushige & Heggiell200ol) . We will not 
discuss the details of the theory here and instead proceed directly 
to the cons e quenc es for the escape rate. This was considered by 
iBaumgardtl ( 1200 ih and he has shown that the relevant time-scale 
for escape is in fact a combination of r^h and Ta 



Trh 



I.e. 



(Al) 
(A2) 

(A3) 



with X ~ 3/4. The proportionality, i.e. the value of the refer- 
ence number A'^i , may be determined by matching the time scales 
following from the theory to the results of A^-body simulations. 
We again define the dimensionless escape-rate ^ as the fraction of 
stars lost per relaxation time such that an additional (small) A'^- 
dependence needs to be included; thus instead of equation il lb we 
have 

l-x 



iV 



3 N_ 

5^ VT-crl 7 \Nl 



The escape rate corrected for the tidal density then becomes 



iVr^r 



71C 



{rh/rj)j 



3/2 



(A4) 



(A5) 



in place of equation BOb . Consequently the A^-dependence also en- 
ters in the expression for x through the constraint that ( is constant 
(equation [lOj 



'^cr'^rh 3 ^ / Tcr 

X= = ^ 

Tcr ^ \ Tcrl 

in place of equation l ll2t . 
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(A6) 



Al Motion in the Tcr — N plane 

The motion in the Tcr — A'^ plane can be found in a similar way as in 
Section l2!4l With A'^ and fcr defined as in equations l lA4b and l lA6b 
we then find 



drcr 



11. 

X Tcr 



iV2 



5 TcrAfl-- - TcrlA^i 
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Because the variables Tcr and N can not be separated, a variable 
substitution u = TcrN^~^ needs to be performed and then the 
variables u and A'^ can be separated and the integration can be per- 
formed. The result is 



Tcr — 



Tcrl 



1 - 



Ao 



-1 5A/2 



(AS) 



Here A = 1 -\- (2/5) (1 — x) and we again used Tcro = 0. When 
X = 1, then A — 1 and equation l lASb reduces to the result of 
Section |2] (equation |19b . 

The continuous growth of Tcr with decreasing A^ in the tidal 
regime shows that clusters are not evolving homologically any 
more in the evaporation dominated regime. This is different from 
the model of Henon where Tcr becomes a constant related to the 
tidal density. The time-dependence also changes and this is dis- 
cussed inlA2l 



A2 Time-dependence 

The time-dependence of the evolution can be found as before from 
an integration over TV (eguation lASb 



N{t) = No(l 

TcvO 



1/^ 



with TcvO now defined as 

(rh/r,,)?/^7V^-- J 



TovO 



Tcr An 



(A9) 



(AlO) 



71C X 

The time-dependent Tcr(t) follows from equation l lASb and the 
above. We can then also get rh(t) and Trh(f), but this will not be 
done here. 



APPENDIX B: TRACKS AND ISOCHRONES IN UNITS OF 

M, ph AND Ra 

In this appendix we present the results of Section|2]in terms of M, 
Ph and pj. We return to the case x = 1. 



Bl Tracks 

The evolution of mass is 

M{Mo,t) ^ Mo -\M\t (Bl) 
If we use m = 0.5 then from equations J16b and J20b we find 



\M\ ~ 35.4C 



_Gpj_ 
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(B2) 
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Together with equations M2\ and l l25t we then have for the density 



Ph(pj,Mo,t) = i (^^^ pj(l- 



. \M\t 
Mo 



5/2 



Ph(Mo,t) ~ i ) , Mo » |M|t, 



Ph(pj,Mo,t) 
or in terms of radius 
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(125Ct)'''^ , Mo > \M\t, (B7) 
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Mo ~ \M\t. 



(B8) 



In the last step we used 14 = 220 kms^^. If we use the constants 
as described in Section|2l (rh/?'j)i = 0.145, C, = 0.2 and G ~ 
4.5 X 10"^ pc^ Mg^ Myr~^ then from equation l lB2l l we obtain 



\M\ ~ 20 Mq Myr"^ ( — 
\kpc 



(B17) 



The tracks and isochrones for ph are shown in Section [3] using 
this isothermal halo approximation. For the isochrones an age of 
13 Gyr was used. 



APPENDIX C: TRACKS AND ISOCHRONES IN UNITS OF 
M, Ph AND Rg, INCLUDING THE ESCAPE TIME EFFECT 

If we include the effect of the escape time (Appendix|All we get for 
the mass evolution 



M(Mo,t)= (Mo"-M,t 
where 



l/x 



M, ~ 35.4C 



^ \ 1/2 
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Here Afi = rfiNi is a constant reference mass that has a similar 
role as A^'i introduced in Appendix [A] and we again adopt rh — 
0.5Mq. Note that Mx is not anymore the time-derivative of M. 



B2 Isochrones 

The isochrones are easily found from the tracks and using A/q = 

M + \M\t 



Ph(pj,M,f) = i (^^) pj(l 



1 + 
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Ph(pj,M) ~ i (^^) pj, M~\M\t, 
or in terms of radius 
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B3 In the isothermal halo approximation 

In an isothermal halo with constant circular velocity Vc the density 
within the Jacobi radius depends on Rq as 



P.i = 



3 

27rG Rl ' 
5.376 f^V 



(B15) 
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CI Tracks 

Using equation l |16t for Tcri, we find that the tracks for general x 
are 

Ph(pj,Mo,t) = -(^-j^p,(^— j h- — 
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(125C^)'''^A^ > lAflt, (B13) rh(pj,Afo,t) 



or in terms of radius 
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C2 Isochrones 



The isochrones for variable x can be found from the tracks by in- 
sertion of an expression for Mo{M, Mx,t) (eguation lClb 
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or in terms of radius 

rh(pj,M,t) = 
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C3 Implications for the comparison to Milky Way globular 
clusters 

The main change that follows from the introduction of a:: 7^ 1 oc- 
curs in the evaporation dominated phase where we do not have a 
constant rh/rj anymore, but instead 

Ph/pj oc M^-"", (C15) 
or 

rh/rj oc M^*"-'^/^ (C16) 

For s = 3/4 and the isothermal halo approximation these scaling 
relations are equivalent to 



Ph oc M'/^Ro\ (C17) 
or 

rh oc M^^^R^J^. (C18) 

So at a given Rq this relation implies that rh oc M^''®, de- 
pending on the exact value of aQ- If we would have taken the 
Coulomb logarithm into account, the index of 1 /6 would be slightly 
smaller. So the particular scaling of Tosc oc t^^'^ has as a result that 
in the regime where mass loss dominates the cluster half-mass ra- 
dius is (nearly) independent of A/ and is determined mainly by the 
tidal field strength. In the expansion dominated phase we still find 
Tcr oc , or ph oc M'^, at a given age. 



^ A relation oc M" is found for x = 1/2. 



